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We discuss opportunities that may arise from subjecting high-multiplicity events in relativistic 
heavy ion collisions to an analysis similar to the one used in cosmology for the study of fluctuations 
of the Cosmic Microwave Background (CMB). To this end, we discuss examples of how pertinent 
features of heavy ion collisions including global characteristics, signatures of collective flow and 
event- wise fluctuations are visually represented in a Mollweide projection commonly used in CMB 
analysis, and how they are statistically analyzed in an expansion over spherical harmonic functions. 
If applied to the characterization of purely azimuthal dependent phenomena such as collective flow, 
the expansion coefficients of spherical harmonics are seen to contain redundancies compared to the 
set of harmonic flow coefficients commonly used in heavy ion collisions. Our exploratory study 
indicates, however, that these redundancies may offer novel opportunities for a detailed characteri- 
zation of those event- wise fluctuations that remain after subtraction of the dominant collective now 
signatures. By construction, the proposed approach allows also for the characterization of more 
complex collective phenomena like higher-order flow and other sources of fluctuations, and it may 
be extended to the characterization of phenomena of non-collective origin such as jets. 



I. INTRODUCTION 

The bulk of low-momentum particles produced in 
heavy ion collisions at the Relativistic Heavy Ion Collider 
(RHIC) and at the Large Hadronic Collider (LHC) ap- 
pears to emerge from a common flow field that is largely 
driven by the density gradients produced in the earliest 
stage of the collision [l|-[4[ . This fluid dynamic paradigm 
of heavy ion collisions is supporte d by a large set of data 
from the LHC and RHIC [IlOf which includes the 
dependence of particle production on the azimuthal an- 
gle relative to the reaction plane and its dependence on 
transverse momentum, pseudorapidity 77, particle species 
and centrality (impact parameter) of the collision. Com- 
plementary support comes from recent indications that 
also event-by-event fluctuations of the energy density dis- 
tribution produced in the earliest stage of a heavy ion 
collision, are propagated fluid dynamically. In partic- 
ular, the measurement of non- vanishing odd harmonic 
coefficients, v n , in the azimuthal dependence of particle 
distributions is an unambiguous signal for the presence 
of sizable fluctuations in the initial stage of a heavy ion 
collision [HH^. 

These findings point to a remarkable set of common- 
alities between the physics of the "Big Bang" determin- 
ing the evolution of the Universe, and the physics of the 
"Little Bangs" formed in heavy ion collisions, that cre- 
ates hot and dense matter under conditions akin to those 
prevailing in the early universe around the first microsec- 
ond, where the transition from the Quark Gluon Plasma 
(QGP) to a confined hadron gas (HG) occurred. 

In both cases, the physical system is viewed at an ini- 



tial time as exhibiting a phase space distribution with a 
high degree of symmetry, overlaid with distributions of 
localized fluctuations. In both cases, these initial fluctu- 
ations are propagated fluid dynamically and are experi- 
mentally accessible, in the first case as fluctuations of the 
(photon) Cosmic Microwave Background (CMB) [2ll-[27j 
and, in the second case, as fluctuations and character- 
istic variations in the density of particles produced in 
very energetic heavy ion collisions |28l |29|. Also, in both 
cases, the decoupling of the different particle species from 
the common fluid dynamical system provides important 
possibilities for experimental verification of the collec- 
tive dynamics. In the case of Big Bang physics, it sets 
the time scale for the decoupling of photons and de- 
termines the primordial abundances of light nuclei. In 
the case of heavy ion physics, it sets, i.a., the hierarchy 
of temperatures for chemical and kinematic freeze-out, 
and determines the relative abundances of hadronic res- 
onances. Furthermore, the study of fluctuations is moti- 
vated mainly by the idea that details of the fluid dynam- 
ical evolution of a system make it possible to constrain 
its material properties. In the case of Big Bang physics, 
this has allowed to constrain the composition of the Uni- 
verse in terms of visible matter, dark matter and dark 
energy [23|. In the case of heavy ion collisions, this has 
allowed to establish that the produced matter exhibits 
the properties of an almost perfect fluid with a shear 
viscosity to entropy density ratio that is close to mini- 
mal UfllailaS. 

The question arises whether these qualitative analo- 
gies between the physics of the Early Universe and the 
physics of heavy ion collisions can provide conceptual 
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or technical advances in either field. Can heavy ion 
physics learn from the techniques developed to analyze 
Big Bang physics, or vice versa, can it contribute to ad- 
vances in modern cosmology? Synergies between both 
fields may arise on different levels. For instance, cos- 
mology has studied in significant detail the question of 
how fluid dynamic perturbations are propagated in an ex- 
panding system, and there is at least some understanding 
of how non-linearities and turbulent phenomena relevant 
for structure formation build up, and which experimen- 
tal measures may give access to them. In the context of 
heavy ion collisions, first analyses of the propagation of 
fluid dynamic perturbations have been undertaken with 
analytical and numerical techniques (32l-[37j , but these 
are arguably at the very beginning, and further explo- 
rations in close analogy to developments in cosmology 
may be beneficial. 

Similarities between the Big Bang and the Little Bangs 
may also lead to synergies on the level of the techniques 
for the analysis of data on fluctuations. The main pur- 
pose of the present paper is to initiate such an effort. To 
this end, we shall use here a standard tool [38| for ana- 
lyzing the CMB anisotropy, and apply it to the study of 
sets of simulated particle densities from ultra-relativistic 
heavy ion collisions. 

We note that studies of the physical system produced 
in the Big Bang and the systems produced in heavy ion 
collisions share not only remarkable commonalities but 
have also significant differences. For an analysis of fluctu- 
ations and their separation from other sources of signals, 
some of the obvious differences (such as the vastly differ- 
ent scales involved in both problems, the different funda- 
mental forces that govern the dynamics, and the differ- 
ence in formulating the dynamics in general relativity or 
special relativity) may be less important. Potentially rel- 
evant, however, is the difference between a study of one 
'big' event (even if consisting of a large but finite sample 
of statistically independent patches) and the study of an 
- in principle arbitrarily abundant - sample of mesoscopic 
events. 

The paper is organized as follows. In section 2 we 
introduce the spherical harmonic approach that under- 
lies CMB data analysis. In sections 3 and 4, we present 
and discuss an application of this approach to simulated 
heavy ion data, and in section 5 we apply the method to 
the particular case of elliptic flow, a standard measure- 
ment in heavy ion collisions. In section 6 we summarize 
our findings and we present a brief outlook. Finally, in 
Appendix A we have presented mathematical details of 
the method and generalization of the model with con- 
stant flow amplitude to more complex one, when the am- 
plitude of the flow depends on pseudorapidity . 



II. FORMULATION OF THE PROBLEM 

In cosmology, the CMB anisotropy is analyzed by 
means of two-dimensional maps that essentially cover 



the full sky (like those from the COBE, WMAP and 
PLANCK experiments [2~0-[23|), or only patches of the 
sky (Boomerang, MAXIMA, CBI, ACT etc. [HQ). 
Analysis of such maps reveal the small temperature fluc- 
tuations and anisotropics that the early Universe has left 
on top of an otherwise smooth and constant background. 
By expressing the measured signal in terms of an expan- 
sion in spherical harmonics, useful quantities such as the 
power spectrum, alignment between different multipoles, 
statistical anisotropy and Gaussianity of the CMB can be 
constructed and compared with theory even from a single 
map (realization) of the CMB sky [21]. In the analysis of 
CMB data, the methods focus on the separation of the 
different components of the observed signal from that of 
the primordial CMB. The by-product of this separation 
yields additional maps for different kinds of foreground 
effects (e.g. synchrotron, free-free and dust emission) and 
the point-sources catalog (39| . 

The analysis of heavy ion data from nuclear collisions 
at the LHC and the RHIC faces similar challenges. On 
one hand, several independent statistical methods are in 
use |40|, |41| to analyze two- and multi-particle correla- 
tions in terms of flow coefficient v n , taken to be sensi- 
tive to the fluid dynamical properties of the high den- 
sity partonic system in its early stages of expansion. On 
the other hand, the same particle correlations are known 
to be sensitive to other important dynamical features of 
heavy ion collisions, such as jets or resonance decays. 
Typical tasks in heavy ion collisions are then to either 
clean the flow signal as far as possible from confounding 
dynamical features, or to analyze a different class of in- 
teresting physical phenomena (such as the 'point-sources' 
produced by jets) without biasing its analysis by under- 
lying collective effects. The question arises, what one 
could learn for this class of problems if the investigation 
of morphological features commonly done in the CMB 
data analysis was applied directly to heavy ion collisions? 

Each heavy ion collision results in a distribution f(r], <j)) 
of particles as a function of azimuthal angle <\> and pseu- 
dorapidity 77. Pseudorapidity is related to the polar an- 
gle of the produced particle with respect to the beam 
direction via r] = — ln(tan(|)), and (j) characterizes the 
angle in the transverse plane orthogonal to the beam. We 
simulate such distributions with HIJING [42| (Heavy-Ion 
Jet INteraction Generator), which is an event genera- 
tor that reproduces a number of main features of heavy 
ion physics in the energy range of RHIC and LHC. It is 
build on PYTHIA [43] routines for hard interactions and 
JETSET [44] routines for string fragmentation. The ini- 
tial geometry and matter distribution is determined via 
a Glauber model [45|, l46j |. Various phenomenologically 
relevant nuclear effects are also implemented, including 
e.g. the nuclear modification of parton distribution func- 
tions, or parametrization of the expected parton energy 
loss. However, HIJING does not directly model a collec- 
tive dynamics resulting in flow. As this is relevant for our 
study, we include flow effects a posteriori by modulating 
the azimuthal distribution of particles from HIJING. 
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Particle production in heavy ion collisions is known to 
display a regular ^-dependence that peaks around 77 = 
and exhibits a forward-backwards symmetry in the case 
of collisions between identical nuclei. Fig. [T] shows re- 
sults for a semi-peripheral Pb-Pb collision at the LHC 
with a total of 12316 particles produced over more than 
16 units of pseudo-rapidity. The pseudo-rapidity distri- 
bution of this randomly chosen simulated single event 
shows visible but relatively small fluctuations around a 
smooth event- averaged distribution. In preparation of 
the following analysis, we have replotted this distribu- 
tion as a function of the polar angle 6. We see that 
the ^-distribution is approximately flat over an extended 
which corresponds to a sizeable win- 
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dow around mid-rapidity (—0.88 < 77 < 0.88). A pre- 




FIG. 1: Top: Pseudorapidity density of charged particles sim- 
ulated by HIJING for a single semi-central collision between 
two Pb nuclei at current LHC energies (^/s nn = 2.76TeV). 
Bottom: The same distribution as a function of the polar 
angle relative to the beam direction. 



sentation in the variable thus highlights the region of 
approximately one unit around mid-rapidity for which 
most differential data are available. 

As remarked already, HIJING does not generate realis- 
tic flow signals. For our study, we superimpose such flow 
signals a posteriori on HIJING events by redistributing 
the generated particles such that event- averaged ensem- 
bles of collisions show a second order harmonic oscillation 
of the azimuthal angle <p with some prescribed magnitude 
V2 . This procedure results in events that carry azimuthal 
oscillations supplemented by event-by-event fluctuations. 
Fig.[2]shows a typical event realization for an average flow 
amplitude V2 =0.10. We caution the reader that experi- 
mental data on flow reveal characteristic dependencies of 
V2 on transverse momentum, rapidity and particle iden- 
tity that are not included in this ad hoc simulation. The 
simplified procedure adopted here will not account for all 
phenomenologically relevant properties of particle flow, 
but it will be sufficient for the purpose of our study. 



III. SPHERICAL HARMONIC 
DECOMPOSITION AND ANALYSIS OF THE 
SYMMETRIES OF THE MODEL 



The harmonic analysis of particle production in heavy 
ion collisions focus typically on the Fourier decomposi- 
tion of the ^-dependence (see Fig. [2j) over a narrow re- 
gion in pseudorapidity around 77 = 0. Here, we contrast 
this procedure with the one employed in the CMB data 
analysis that starts from an expansion of the entire two- 
dimensional map /(#, (j)) in terms of spherical harmonics 
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FIG. 2: Density of charged particles as a function of the az- 
imuthal angle <j) in the plane transverse to the beam direction, 
simulated by HIJING and modified by an elliptic flow with 
v 2 0.10. 
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where a^ m = |aj ?m | exp(— im<&i, m ) are the coefficients of 
decomposition with amplitudes \a^ m \ and phases ^/, m 
for each component Z,m, Pj m (cos0) are the associated 
Legendre polynomials, and 
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is the total power spectrum. 

As it is seen from Eq. (pQ) , the constant term with I = 
which corresponds to the average density of the parti- 
cles on the sphere is not included in our analysis. Thus, 
(j)) represents the fluctuations of the number density 
of particles with respect to the mean value and it can 
take positive or negative values. The parameter l ma x in 
Eq. (pQ) characterizes the angular resolution of the map. 
In principle, infinitely precise pointing resolution requires 
(Imax —> oo ) which is computationally prohibitive. In the 
CMB analysis, one choses the parameter l max accord- 
ing to the angular resolution of the experimental data 
which is dictated by the segmentation of the detector. 
For instance, l max = 100 corresponds to an angular res- 
olution A = 7r /lmax — 1.8°. In heavy ion physics, other 
physics considerations such as the typical angular sepa- 
ration of nearest charged tracks in the Arj and A(f> may 
motivate the choice of Z max . In the following, we use a 
resolution lmax = 100 that is sufficiently high to reveal 
fluctuations on small scale. We work for illustrative pur- 
poses with the standard GLESP [38| pixelization where 
lmax = 100 corresponds to a number of pixels in the 0- 
direction Nq = 2l ma x + 1 and in the azimuthal direction 
= 2N e . 



A. Mollweide projection of heavy ion events 

To familiarize ourselves with the main features of a 
CMB-like harmonic analysis of heavy ion collisions, we 
start from a 'baseline' distribution /(#, </>) of a heavy ion 
collision in which features of global collective dynamics 
are absent. To this end, we consider as baseline a HI- 
JING event without superimposed flow signal (v2 = 0). 
Fig. [3] shows the corresponding fluctuations of the distri- 
bution /(#, (j)) in the so-called Mollweide projection that 
is heavily used in CMB analysis. In this representation 
the polar angle 6 runs vertically from the 'north pole' to 
the 'south pole' and the azimuthal angle (j) runs along the 
'equator' from the center to left . Let us pause to discuss 
in more detail the information shown in this map. 



We have seen already in the context of Fig. [T] that 
particle distributions are approximately flat in a wide 
window of the polar angle j < 6 < Accordingly, 
the particle density and fluctuations in remain approx- 
imately unchanged in a broad band around the equator 
of the Mollweide projection (see upper part of Fig. [3j). 
Larger particle densities are found at the poles reflect- 
ing the large number of particles emitted at small angles 
with respect to the beam direction (see bottom panel of 
Fig. [1]). This enhanced particle yield at the poles is a 
kinematically trivial but potentially confounding factor 
for studies that aim at establishing event-by-event global 
or local signals on top of an event-average background 
that varies significantly in r] or 6. Thus, irrespective of 
the choice of coordinates, the question arises to what 
extend potentially interesting information about fluctu- 
ations and azimuthal asymmetries can be disentangled 
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FIG. 3: (Color online) Top: Mollweide projection of a single 
HIJING event without flow. The azimuthal angle (f) runs along 
the 'equator', the polar angle runs from pole to pole. The 
color coding tracks deviations relative to the average level of 
counts. Bottom: The m — mode has been removed from 
the map. The angular resolution of the maps corresponds to 
lmax = 100 (see Eq. (pQ). The multiplicity of the event is 
17000. 
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from the strongly varying background in rj or 6. The 
approximately uniform distribution of fluctuation in the 
lower panel of Fig. [3] demonstrates that removing from 
the expansion of /(#, <j)) the m = mode of the har- 
monic decomposition goes a long way towards achieving 
this aim. We note that the m = spherical harmonic 
components do not depend on <j) but only on 9. There- 
fore, removing this component will remove the dominant 
dependence on the average pseudorapidity particle distri- 
bution without modifying any flow signal. In general, a 
dependence of particle distributions on pseudo-rapidity 
will remain after subtraction of the dominant m = 
component, of course. In appendix A, we explain how an 
arbitrary ^-dependence can be dealt with in the present 
formalism. For the sake of a particularly transparent pre- 
sentation, we focus in the main text on the simpler case 
in which the ^-dependence can be treated satisfactorily 
by subtracting the m = mode. 

To further quantify this separation of background from 
fluctuations, we discuss now in more detail the spherical 
harmonic analysis of 0). Using Eq. (pQ) and, for in- 
stance, the GLESP transform method, we can decompose 
this distribution into its multipole components (see Ap- 
pendix A for details). The result of this operation is the 
set of amplitude coefficients a^ m . With a suitably high 
choice of l m ax all morphological features of the initial 
angular distribution of counts are preserved. For the de- 
composition we have used I max = 100 in accordance with 
the original binning of the fluctuations shown in Fig. [3j 
The corresponding power spectrum C(l) is shown in the 
top panel of Fig. 2J 

In the case of a perfect reflection symmetry of particle 
production around mid-rapidity, — § <H> — (0 — ^), the 
odd components C s (l) of the total power spectrum van- 
ish. Event-by-event fluctuations, however, will break this 
reflection symmetry. As a consequence of such fluctua- 
tions, the top panel of Fig. 3] shows very small but non- 
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FIG. 4: (Color online) Top panel: power spectra for the map 
from Fig. [3] The black line is the total power spectrum 



vanishing values C s (I) G [10 5 ; 10 4 ] for essentially all C s (l), the blue dots are for the m = mode only and the red 



odd integers I. In contrast, the even modes of C s (I) are 
for small even integers I up to a factor 1000 larger since 
they are dominated by the strong dependence of parti- 
cle production on pseudorapidity. Removing the dom- 
inant m = mode from this power spectrum, we see 
that all even and odd modes in C s '(I) have now approx- 
imately equal signal strength, and that the odd modes 
do not change in strength. This quantifies the extent 
to which removing the m = mode removes the trivial 
background dependence without affecting information on 
fluctuations. We conclude that this procedure allows one 
to represent a 'baseline' event, in which fluctuations are 
superimposed on top of a varying background, in a way 
in which the varying background is efficiently removed 
and fluctuations dominate the visual representation (see 
Fig. [3j) and the statistical information (see Fig. [4]). We 
shall discuss in the following section how this changes if 
the event contains, for example, a collective flow compo- 
nent. 

We finally mention yet another way of checking that 



line is for the \m\ > 1 modes. Bottom panel: the normalized 
distribution function H(s) for versus amplitude of fluctuations 
s for the total signal (the black line) and for the signal with 
| to | > 1 (the red line). 



removal of the m = mode isolates information about 
the nature of fluctuations in the event. To this end, we 
recall that m = (— l) m a^_ m , and we write the power 
spectrum of Eq. ([2]) in the form 



C s (l) = C(l)+D(l) 



\ a l> 



2/ + 1 



21 



aw. 



(3) 

The first part of Eq. ([3]) , C(Z), corresponds to the 
most symmetric m = azimuthal modes, while the sec- 
ond part, D(l), shows the power spectrum of fluctuation 
above the \m\ = threshold. It would be worth to note, 
that both these components have significantly different 
statistical properties, which can be characterized in terms 
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of the probability distribution functions. The bottom 
panel of Fig. 3] shows the difference in the statistical 
properties of the signal when all m modes are included 
(black) and the case \m\ > 1 (red). Normalized to the 
total number of pixels ( N to t — ^ ma x)^ the probability 
distribution function, H(s), is defined as the number of 
pixels with corresponding amplitude of fluctuation within 
an interval [s-Ss, s + 5s]. Here, 6s = (s max - Smi n ) /Nun, 
s max and are the absolute maxima and minima of 
the signal in the map with m > 1 (red), and A/^ n = 200 is 
a number of intervals ( bins). As it is seen from Fig. |4]the 
probability distribution function for signal with removed 
the m = mode, is quite close to a Gaussian distribu- 
tion, while the one including all m modes (black) is far 
from a symmetric Gaussian distribution with respect to 

5 = 0. 

To the extent to which the m = 0-mode dominates 
the distribution, any distortion of the morphology of the 
map will provide relatively small corrections to the total 
power spectrum C 5 , but it may result in strong modifi- 
cations of the power spectrum D(l) (see Eq. (j3j)), and of 
the corresponding multipole coefficients a^ m . Note that 
any features of the a^ m and of the corresponding power 
spectrum can be detected only if they exceed the level of 
statistical event-by-event fluctuations. This is why the 
above mentioned 'Gaussianisation' trend of the probabil- 
ity distribution function for the sub-dominant component 
power spectrum, D(l) , reflects the level of detectability 
of particular features with small amplitudes above the 
statistical level. 



IV. SINGLE EVENT v n -FLOW 

Due to Lorentz contraction the colliding ultra- 
relativistic nuclei appear pancake-like just before a heavy 
ion collision in the laboratory frame of reference. For 
non-central collisions, the overlap between the two nu- 
clei has an approximately almond-shaped cross section. 
This initial and anisotropic collision geometry results in 
large pressure gradient differences along the principal 
axes of the distribution and gives rise to a fluid dynam- 
ical flow characterized by fundamental properties of the 
fluid like the reinteraction time and residual interaction 
between the constituents. This flow, in turn, leads to an 
anisotropic distribution of particle momenta and num- 
bers in the transverse plane which may be quantified as 
a function of the azimuthal angle, (j) (lol l4lj. 

Recent experimental results from both RHIC and the 
LHC suggest that the particle flow is established at the 
quark-gluon level over a characteristic time scale of about 
1 — 2/m/c and that the flow is quite sensitive to detailed 
features of the system, such as the viscosity [30|, |3l| . Re- 
cently, the presence of higher order flow moments was 
understood as resulting from fluctuations around the al- 
mond shape in the initial collision geometry fl^-[2Q| . 

In the present study, the angular distribution of parti- 
cles produced in each heavy ion collision will be viewed 



as a 'stochastic' map 



s(M) = /(M) i 



2 v n cos[n( 



*»)] 



(4) 



This expression maps the random distribution of parti- 
cles without flow, f(6, <p) = d 2 N /d(j)dr]\ Vn= o onto a distri- 
bution S(9, <j)) = d 2 N/d(j)drj that shows azimuthal depen- 
dencies with harmonic flow amplitudes v n in azimuthal 
orientations \I/ n . The model distributions considered in 
the following were obtained in line with this expression 
by modulating, a posteriori, an azimuthal uniform distri- 
bution f(6, (j)) simulated in HIJING or a simple random 
generator with prescribed flow amplitudes v n with angu- 
lar orientations \I/ n 1 . The case n = 2 corresponds to 
so called elliptic flow, to which we restrict ourselves in 
the present study. For simplicity, the model used in this 
section assigns the same flow value to all regions of n and 
transverse momentum p t , and to all particle species (e.g. 
mesons and baryons). This could be extended easily to 
account for dependencies in these parameters, or to treat 
the flow coefficients as random functions with statistical 
properties that differ from those of f(9,(j)). The argu- 
ments in the present study will not require such a more 
detailed modeling. 

A typical experimental task is to extract from a given 
particle distribution d 2 N/dcj)dr] the flow coefficients v n , 
defined as [47j 



M = ((e 



» 



(5) 



Here ((..)) denotes the average over all particles for each 
heavy ion collision and over the entire statistical ensemble 
of events. The following discussion will not involve this 
later step. Rather, we shall discuss how different flow 
coefficients v n and event planes \I/ n can be identified in 
a CMB-like harmonic analysis of single heavy ion events. 
To this end, we decompose Eq.(J2|) in spherical harmonics 



.(6) 



Here, are the coefficients of the spheri- 

cal harmonic decomposition for S(9,(j)), and 
Q,mTn a>i, m g(hm,Tn), where g(l,m^n) 

2nN hrnTn N hm J 1 _ 1 dxPj n ^ n (x)P{ n (x). This type of 
equation is also encountered in the CMB data analy- 
sis [48j for cases where the statistical isotropy of the 
signal is broken by regular modulations. We now discuss 
this equation in more detail. 



1 In general, the event plane angles ^f n are correlated to the az- 
imuthal orientation of the reaction plane, that is denned as 
the plane spanned by the impact parameter vector and the beam 
direction. However, while and the ^ n 's are correlated, they 
are in general not identical and different harmonic modes can be 
correlated differently. Therefore, there will be in general an event 
plane angle ^ n associated with each harmonic n, that may be 
determined in a more comprehensive analysis in terms of higher 
order harmonics, and different ^ n 's will not coincide in general. 
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FIG. 5: (Color online) Schematic representation of the contri- 
bution of V2 (elliptic) and V4 flow to the various m components 
of the &z = 4 ?m -coefficient. The ^-modulation will change the 
amplitude and the phase of the 64,2 coefficient without con- 
tributing to the other components. The va modulation will 
only change the 64,4 component. 



FIG. 6: (Color online) The power spectrum for an event with 
elliptic flow V2 — 0.07. The total power C s (l) is shown in 
black, the blue stars correspond to C s (I) — D(l) = C(Z), and 
the D(l) power spectrum is shown in red (see Eq.(3)). 



A. v n - modulation with n — even 

The iterative solution of Eq. (j6j) has an especially sim- 
ple form for even-n, n = 2k, k = 1,2... In this case, 
the dominant component of the a\ ^-coefficients is asso- 
ciated with m = modes, and I = even. As we have 
pointed out in Section 3, the strongest component of the 
f(0, (j)) signal is the m = component (blue dots in Fig. 

This implies that, in Eq. ©, the last term in the 
brackets a^ m _ n is maximal, if m. 

For m = 0, Eq. ([6]) reads 

bi,o = a,, + J2 u » (d,»e~* n * B + c (> _„e in *") 

n 

= a z>0 + 2 y^Vn\ci, n \ cos(n(<fo >n - ^ n )) , (7) 



coupling, but this effect is significantly smaller (~ V2), in 
comparison with the discussed 64,2 and 64,4 components. 
These features of Eq. ([6]) motivate to use 



K,n = cin.n + v n b n $g (n)e m * n , 



(8) 



as a basis for separating the contributions from differ- 
ent v n , where g(ri) = 2itN n $N nin dxP®(x)P™(x) (see 
Appendix A). 

For a sufficiently large flow signal, such that |a n5n | <C 
v n \bn,o\g(n)i the iterative solutions of Eq. (j8j) for the am- 
plitude of the flow v n and the event plane \I/ n for a single 
event are then given by 



K,n\ .j. , 

9{n)\on,o\ 



(9) 



where we have used the relation a/ 5 _ m = (— l) m a z * m . For 
the scenario studied here, we can neglect in Eq. ([7]) the 
second term (ex v n \c^ n \) relative to the first one, due 
to the inequality \a^o\ ^> \ c l,n\ which follows from the 
dominance of the m = modes. We consequently obtain 
the following approximate solution: b^o ~ a^o. 

The coefficients bi jTn vanish for I < \m\. As a conse- 
quence, in Eq. (j6j), for instance, the coefficient bi y7n with 
m — 2 and I = 2 is connected to the V2 component only, 
and does not receive contribution from higher flow mo- 
mentum with n > 2. For I = 4 there exists only two 
terms in b^ m which are in resonance with the 64^0 mode: 
^4,2 through the V2 term in Eq. ([6]), and 64,4 through 
V4. We illustrate these relations further in Fig. [51 We 
note that the elliptic flow will also provide a 64,3 o 64,1 



B. v n - modulation with n — odd 

The method presented above clearly illustrates how 
spherical harmonics with I = even give access to am- 
plitudes of the v n -&ow coefficients with n = 2,4,6... We 
now focus on the odd harmonics v n , n = 1,3,5,... that 
are known to provide information about the fluid dy- 
namic response to fluctuations present in the initial stage 
of heavy ion collisions. To this end, we study the &z >m - 
basis of Eq. ((6]) for the case I = even and n = odd. We 
start from the most powerful I = 2 and I = 4 compo- 
nents of b^ m under consideration (see Fig. [4]). For the 
quadrupole component, I = 2, in the presence of even 
and odd flow harmonics, the "resonant" component 62, ™ 
for v\ is 62,1' 



where the phase of the event plane corresponds to 
^i-flow. In Eq. (fTU|) the major contribution to 62,1 is 
given by the first term, proportional to the most powerful 
component <22,o — &2,o? and 



&2,i - vih,ogie 



Vl 



gi\h,o I 



^2,1 



(11) 



where 02, 1 is the phase of the 62,1-component, and g\ = 
47r J Q (ixP 2 3 (x)P 2 1 (x). It is obvious that for any odd n the 
corresponding solution for the amplitude and the phase 
is given by 



+ l,n| 



g(n+l)\b 



'n+l,0| 



£>n+l,n • 



(12) 



where g(n + 1) = 4tt ^P° +1 (x)P^ +1 (x). Similar to 
n = even case, the phase of the event plane for n = odd 
can be reconstructed also from <p n ^ n phase as follows: 
n^ n = n>n . 

Note that all methods of reconstructing the orientation 
of the reaction plane (s) in heavy ion collisions distinguish 
between the 'true' reaction plane Sfr n and the orientation 
of the event plane reconstructed from experimental data 
. Uncertainties in extracting the true parameter 
from a statistically finite amount of data are character- 
ized by quoting the finite resolution. That means that, in 
respect to the statistical ensemble of realizations, v n and 
\I/ n should be treated as random variables with probabil- 
ity density functions P Vn and P^ n . Performing the same 
estimation as in Eq. (J9j) and Eq. ([T2]) for each event 
we can determine the number of realizations with v n and 
\I/ n within a given interval of uncertainty, and in fact, we 
can obtain the probability distribution functions P Vn and 
P^ n for a given ensemble of realizations. Then, one can 
define the first moments of the corresponding variables 
(v n ) and (^n), for a statistically finite amount of data. 
We will illustrate this approach in the next section. 



V. ELLIPTIC FLOW: THE CASE n = 2 

Here, we further illustrate how the symmetries of the 
spherical harmonic representation of the particle distri- 
bution can be used to extract the elliptic flow signal fol- 
lowing Eq. (|9]). As discussed in Section 3 and shown in 
Fig. O this signal is a modulation of the random back- 
ground particle distribution by a factor ex v n cos(2(/> — 
In Fig. [6] we show the corresponding power spectrum us- 
ing a flow of V2 = 0.07. The power in the components 
D(l), I = 2,4,6,8 is seen to be significantly higher than 
the asymptotically flat and small D(l) ~ const values, 



0.03+0.00- 




FIG. 7: (Color online) From the top left down to the bot- 
tom right: Maps of the 1 622 1 amplitude for input flow values 
of v 2 = 0,0.01,0.05,0.075,0.1,0.2,0.3,0.5, respectively. For 
simplicity = in all maps. A different reaction plane an- 
gle would manifest itself as a phase shift in the 'East- West' 
direction. 



obtained from a random realization without flow (see Fig. 
[4j). It is in this way that the CMB-like harmonic analy- 
sis identifies - on the basis of a single event and without 
recourse to an event sample - the presence of a signal on 
top of random localized fluctuations. The enhancement 
of the even low-1 multipoles originates from flow contri- 
butions to 62,2, &4,2? &6,2- Expressing these multipoles via 
Eq. (j6]), we find 



D(l) 



+ 1 ^ 

m=l 

v 2 2 C(l)g 2 (n)+D(l), 



21 



J rv 2 bi^g{n)5 



^ |2 



(13) 



where (21 + l)D(l) = £ m=1 \a hm \ 2 . 

We discuss now in further detail the relation between 
the flow signal and the multipole components determined 
in a CMB-like harmonic analysis. In a Mollweide projec- 
tion, one can visualize easily how the 1 62,2 1 component 
changes with a flow signal varying between v 2 = to 
V2 = 0.5, see Fig. [71 In Fig. [5J we demonstrate this de- 
pendence beyond the level of single events by averaging 
over samples of 10 5 events with a specific flow amplitude 
V2. We see that there is a linear relation between v 2 and 
the ratio of the dominant quadrupole components |&2,2 1 




"D 



0.5 



FIG. 8: (Color online) Estimator D = |6 2 ,2 1/^(2) |& 2 ,o| for 
the elliptic flow amplitude obtained for single events plotted 
versus the V2 value used in the simulations. The error bars 
correspond to 68% confidence level. 



and 1 62,0 1- The small error bars shown in Fig. [8] corre- 
spond to a la standard deviation. They illustrate the de- 
gree to which in the present model study fluctuations at 
the single event are small compared to the event- averaged 
mean. An analogous conclusion can be drawn from Fig. [9] 
about the correlation between the true event plane (^r), 
which is known for each simulated event, and the phase 
(/>2,2 determined in the analysis (see Eq. JSJ). The in- 
sert shows the distribution of the difference — 0.502,2 
which is a measure of the event plane resolution obtained 
in this method. The distribution is consistant with a 
Gaussian with mean value \i ~ and standard devia- 
tion a = 0.0254rad. The tight correlation between both 
orientations is due to the fact that particles at all rapidi- 
ties show on average the same azimuthal correlations and 
this contributes to constraining statistical fluctuations in 

02,2- 

We finally indicate how a CMB-like harmonic analy- 
sis can be used to assign to an event a probability that 
the measured harmonic components result as fluctuations 
from a random background, rather than in response to 
a collective flow signal. To assess this question, we have 
generated the distribution of |a2,2 1 values obtained by 
decomposing 10 5 events without flow into spherical har- 
monics. The probability P of assigning a given event as a 
random background for, say V2 =0.01, can be defined as 
the number of realizations without flow above the limit 
of 1 62, 2 1 corresponding to V2 = 0.01. We have found 9230 
out of 10 5 random realizations, and the probability of 
correctly concluding an elliptic flow value of V2 = 0.01 to 
an event with this |&2,2 1 analysis is therefore 1 — P ~ 90%. 
For a five times larger signal V2 = 0.05, the corresponding 
level of detectability is better than 99,999%. We paral- 
lel this analysis for the harmonic component |&2,i| an d 




r 22 



(rad) 



FIG. 9: (Color online) The reaction plane angle (rad) for 
each of 202 HIJING events with flow = 0.07 versus the phase 
02,2 (rad) obtained from the analysis described in Eq. ©. 
The insert shows the distribution of the difference ^^—0.502,2 
(rad). 



have confirmed the fact that the limits for V2 = 0.01 and 
V2 = 0.05 both lie around the average of the background 
distribution. 

The exploratory analysis presented here involves infor- 
mation about the simulated events that is not directly ac- 
cessible experimentally. For instance, experimental data 
do not allow to check directly the correlation between 
1 62, 2 1 an d the 'true flow signal i^', or between ^2,2 and 
the 'true reaction plane as done in Fig. [8] and Fig. [9] 
In this sense, our discussion serves mainly the purpose 
to illustrate that a CMB-like harmonic analysis can be 
applied to individual heavy ion events and that it has the 
sensitivity of associating flow values and event planes to 
single events in a meaningful way. We further note that 
while Figs. [8] and [9] are not directly measurable, closely 
related quantities are. For instance, one could plot the 
ratio of the quadrupole components |&2,2 1 and |&2,o | as 
a function of event centrality to characterize the nature 
of event-by-event fluctuations in an experimentally mea- 
sured event sample. We also note that we have tested the 
sensitivity of the method to a coarser binning of the pri- 
mary particle multiplicity information, as it could result 
from a finite binning of the detectors. We find no signif- 
icant difference in the results when using 100 bins in rj 
and 40 bins in <j>. Also, methods have been developed for 
the CMB analysis to handle partial detector acceptance 
and efficiency etc. A more detailed discussion of these 
aspects will be left to further work. 
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VI. SUMMARY AND CONCLUSIONS 

The present work is a bold exploratory step aimed at 
testing the use of CMB-like analysis tools for the study of 
high- multiplicity heavy ion collisions. It parallels other 
recent efforts, see e.g. [28|, [29|. We started from repre- 
senting simplified simulations of heavy-ion collisions in a 
Mollweide projection commonly used in the CMB anal- 
ysis. We then asked how pertinent features of heavy ion 
collisions, including collective flow and small-scale fluc- 
tuations manifest themselves visually in such a represen- 
tation, and how they are characterized by standard CMB 
tools. 

A CMB-like analysis of fluctuations is based on an ex- 
pansion in terms of spherical harmonics and thus returns 
the two-parameter set of coefficients fy ?m . For the charac- 
terization of collective flow phenomena, that show up as 
purely azimuthal dependencies, the azimuthal harmon- 
ics with flow amplitudes v n provide a complete set of 
functions. Therefore, the larger set of coefficients bi y7n 
must arguably contain redundancies for the characteri- 
zation of collective flow in heavy ion collisions. However, 
event- wise fluctuations in heavy ion collisions are not con- 
fined to the azimuthal (^-direction and may be expected 
to show locally no preferred orientation in the ?7-</>-plane. 
Our study indicates that the larger set 6^ m , while show- 
ing redundancies for the characterization of ^-dependent 
phenomena, may offer novel opportunities towards char- 
acterizing small-scale fluctuations in heavy ion collisions 
and separating them from flow effects. For instance, we 
observed that due to the r] —> — r] reflection symmetry 
of event- averaged identical nucleus- nucleus collisions, ex- 
pansion coefficients with even integers I are inevitably 
dominated by event- wise fluctuations. The set of these 
coefficients may then be used as a baseline on top of 
which collective or global event properties can be estab- 
lished. We have discussed all possible generalizations of 
the method in Appendix A, focusing on ^-dependency 
of the amplitude of the flow .This baseline allowed, for 
instance, to check in Fig. 2] to what extent removal of 
the m = mode subtracts efficiently the effects of a 
rapidity-dependent multiplicity distribution from a fluc- 
tuation analysis. And we have seen in the discussion 
of Fig. [6] how the larger set of coefficients bi jin can help 
visually and statistically to identify flow on top of a fluc- 
tuating background. 

For the detection of flow-like azimuthal modulations, 
the method used here is limited by the level of statistical 
noise v n > N~z^ where N is the multiplicity. Thus, for 
v n ~ 10 ~ 2 the corresponding multiplicity required is in 
order of 10 4 particles, which is close to the parameters 
at the LHC. This is parametrically the same dependence 
as e.g. in a flow analysis in terms of second order cu- 
mulants. Higher order cumulant expansions can improve 
this statistical limit. In this sense, the CMB-like analysis 
of heavy ion collisions shown here provides an alternative 
for characterizing collective flow but does not offer obvi- 
ous parametric advantages. However, it offers arguably 



a very well-suited setting for analyzing those event-wise 
fluctuations that remain once the dominant azimuthal 
modulations have been removed. In analogy to the ques- 
tions asked in cosmology, one can then characterize what 
sets the scale of these remaining fluctuations. For in- 
stance, which physics underlies the fluctuations seen in 
a Mollweide projection of heavy ion collisions after re- 
moval of flow harmonics, such as in Fig.[3j How are these 
fluctuations sensitive to resonance decays, jet-like parti- 
cle correlations or jets? On which scales could critical 
phenomena leave characteristic signatures? Which simi- 
larities or characteristic differences could be expected for 
the corresponding maps based on fluctuations in the dis- 
tribution of electric charge, or strangeness content? We 
intend to follow up these and further questions in the 
analysis of real data, as well as in further model studies. 
Some of these questions are well known in CMB science, 
when we can use , for instance, optimal filter approach to 
amplify the point-like sources , incorporated into diffuse 
background. This method in combination with wavelet 
transform of the signal can potentially detect localized 
features with very small amplitude of heavy ion collisions, 
which can be associated with low amplitude jets. This 
work is in progress and will be published in a separate 
paper. 
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Appendix A: Stochastic equation and general 
description of modulation 

In this appendix we discuss how the method of sec- 
tion IV can be generalized to an analysis that accounts 
for an arbitrary dependence of v n = v n {0) on pseudo- 
rapditity r] = — ln(tan(|)). This can be done by mul- 
tiplying the stochastic equation by functions W(ff) with 
different shapes. The stochastic equation for modula- 
tions of the random signal is given by: 



S(8,<f>) = f(0,< 



l + 2^2v n cos(ra(<£ - *„)) 



]T v n [e in *"s+(6, 4>) + e- in ^s-(0, cj>)] 



(Al) 
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where 

s+(e,<j ) ) = f(e,<j>)e- in t s-(e,4>) = f(e,4>)e in<i >. (A2) 

The distribution S(9,(j)) can have a significant 0- 
dependence. We weigh the left and right hand side of 
Eq. (|Al|) by the window function W(0), 

FS(0,(/>)W(0) = FS(0,(j))W(0) + 

v n FW{0) [e in *"s + (0^)+e- in *"s-(0^)] ,(A3) 

where F is the linear operator of the spherical harmonic 
decomposition 



n27T f'TT 

\.. = d4> 

JO JO 



d0 sin 0Y l * rn (0,< 



(A4) 



In radioastronomy, this is referred to as the apodization 
approach and it serves to amplify some particular range 
of the (0, (j)) map, in order to reduce some very bright 
sources of contamination of the signal under investiga- 
tion. Here, we use it to select bins in for a characteri- 
zation of the dependence of the signal on pseudo-rapidity. 
The coefficients of the spherical harmonic decomposition 
read then 

k,m = Cl,m+Vn / / d0 sin X 

Jo Jo 

k, m = FW(6)S(e,<t>), c l:m = FW(6)f(e,<t>). (A5) 

Using the spherical harmonic decomposition /(#, (j)) = 
Ez',m' a ^m'^',m'(0,0), from Eq(|A5j) we have 

ci, m = FW(0)f(0^)=J2 a i',rn>Gl>™'(n = O), 

V ,m' 

I' ,m' 



(A6) 



I'm' 



where 



G l i;™\n) = dcji [* <w emerge,*) x 

Jo Jo 

ypyiM)^. (A7) 

The coefficients of the spherical harmonic decomposition 
satisfy then 



C b l,m ^ C b 



Lm 



(A8) 



It is technically advantageous to introduce the variable 
x = cos 0. With the help of 



I (21 + 1) (I - m)\ 
Air (I + m)\ 



(A9) 



one can then express the coefficients in Eq. (|A8|) as 



l',m' 
— ^m,n 

ai '>° Gl l,n( n ) + °( a l,o), 

V 

G l ^(n) = 2nNi t0 N Vtn x 

y dxW(x)P^(x)P^ n (x) (A10) 

The case of a ^-independent distribution that we dis- 
cussed predominantly in the main text, is recovered for 
W(x) = 1. In this case, one can use the integrals 



1 dxp?{x)p?{ x ) = m ]) 



21 + 1 ' 

1 dxlP m (x)] 2 _ 2 (l + m)\ 
dx\P x (x)\ - 2l + l(l _ m)l 



to write 



G 2 ; 2 (n = 2) = --j=, G 2 ] = l. 
According to Eq. (|6j), this gives 



V2 



"2,0 



(AH) 



(A12) 



(A13) 



In order to analyze the ^-dependency of the amplitude 
of the flow v n (0), keeping in Eq. (|Al|) the azimuthal mod- 
ulations proportional to cos(n(</> — \£ n )). Then, we can 
decompose the amplitude of the flow through Legendre 
polynomials: 



Vn(0) 



^Vc;p g (cos( 



q 



(A14) 



where: v n = const, and are the coefficients of de- 
composition. Obviously, the decomposition of v n (0) in 
the form of Eq. (|A14p is not unique. One can use, for 
instance the associated Legendre polynomials P{ L (0) 1 in- 
stead of Legendre polynomials. The particular choice of 
the orthogonal functions should reflect the properties of 
the effects of modulations under investigation. 

In the most general case, azimuthal modulations are no 
longer separable from polar ones. In heavy ion physics, 
this is the case e.g. for jet-like particle correlations. In- 
stead of Eq. ([Al|) , the starting point for analyzing such 
modulations would be then the stochastic equation 



S(0, </>) = f(0, <f>)[l + 2V(8, </>)] , (A15) 



where 



(A16) 



I m— — l 
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The coefficients S^ m of the spherical harmonic decom- 
position can be expressed as Si jTn = a^ m + D^ mi where 

oo V oo i" 

A,™ = 2^ E E E (-^"flwwx 

Z'=l m'=-V l"=0 m" = -l" 

/ (2Z" + 1)(2Z' + 1)(2Z + 1) ( I" V l\ 
X V 4^ ^ Oj X 

x ( " ^ ' I (A17) 
I ra" m' — m y ' v 7 

and a/ m = Ff(6,(j)). This equation is well known in 
CMB physics (see for instance (49j| ) . The major differ- 



ence betwen modulation of the signal in form of Eq. (|A16|) 
and v n — const is related to redistribution of the power 
from even multipoles to odd ones, if Vbn+i,™ has non-zero 
components. In the case when only V2n,ra-coefricients 
are non-vanishing, this quadrupole modulation will re- 
distribute the power from the n-th mode to n — l,n + 1 
modes, where I = 2n and change the orientation of the 
quadrupole from very planar ( V2 = const) to non-planar, 
depending on amplitude of 2, 1-component. We shall ex- 
plore in future work to what extent Eqs. ()A16|) and (|A17p 
allow to investigate more complex modulations of the 
particle distribution in high-multiplicity heavy ion col- 
lisions. 
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